---
title: "JOP_repl"
author: "Olga Chyzh"
date: "9/3/2019"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,message = FALSE,warning = FALSE,fig.cap = "...")
```

Libraries:
```{r}

library(dplyr)
library(ggplot2)
library(maps)
library(mapdata)
library(leaflet)
library(mapproj)
library(tidyverse)
library(magrittr)
library(scales)

```

Open the data:
```{r}
rm(list=ls())
trumpilot<-read.csv("soydata.csv", header=TRUE)
fixeffs<-read.csv("fixeffs.csv", header=TRUE)
```


Make the base border of the US map:
```{r}
counties <- map_data("county") #map_data function is under the ggplot2 package. 
states<-map_data("state")
p <- ggplot() +
  geom_polygon(data=counties, aes(x=long, y=lat, group=group), fill="grey85", color="white")

usamap <- p+
  geom_path(data=states, aes(x=long, y=lat, group=group), color="black")+
  theme(axis.text=element_blank(),
              axis.ticks=element_blank(),
              axis.title.x=element_blank(),
              axis.title.y=element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.border = element_blank(),
              panel.background = element_blank())

```

Merge with county border data:
```{r}
counties$subregion<-str_remove(counties$subregion," county")
trumpilot$subregion<-str_remove(trumpilot$subregion," county")
trumpilot$subregion<-str_remove(trumpilot$subregion," parish")
trumpilot$subregion[trumpilot$subregion=="dekalb"]<-"de kalb"
trumpilot$subregion[trumpilot$subregion=="desoto"]<-"de soto"
trumpilot$subregion[trumpilot$subregion=="dewitt"]<-"de witt"
trumpilot$subregion[trumpilot$subregion=="dupage"]<-"du page"
trumpilot$subregion[trumpilot$subregion=="hampton city"]<-"hampton"
trumpilot$subregion[trumpilot$subregion=="lamoure"]<-"la moure"
trumpilot$subregion[trumpilot$subregion=="laporte"]<-"la porte"
trumpilot$subregion[trumpilot$subregion=="lasalle"]<-"la salle"
trumpilot$subregion[trumpilot$subregion=="newport news city"]<-"newport news"
trumpilot$subregion[trumpilot$subregion=="norfolk city"]<-"norfolk"
trumpilot$subregion[trumpilot$subregion=="o'brien"]<-"obrien"
trumpilot$subregion[trumpilot$subregion=="prince george's"]<-"prince georges"
trumpilot$subregion[trumpilot$subregion=="queen ann's"]<-"queen anns"
trumpilot$subregion[trumpilot$subregion=="suffolk city"]<-"suffolk"
trumpilot$subregion[trumpilot$subregion=="virginia beach city"]<-"virginia beach"
trumpilot$subregion[trumpilot$subregion=="yellowstone"]<-"yellowstone national"
trumpilot$subregion<-str_remove(trumpilot$subregion,"\\.")

counties<-counties %>% mutate(subregion=trimws(subregion), region=trimws(region), id=paste(region,subregion, sep="_"))
trumpilot<- trumpilot %>% mutate(subregion=trimws(subregion), region=trimws(region),id=paste(region, subregion, sep="_"))
mydata<-left_join(counties, trumpilot, by=c("region","subregion"))

```

Figure 1: Change in Vote Share Between the 2016 and 2018 Congressional Elections
```{r}

mydata$votech2<-NULL
mydata$votech2[mydata$votech/100<(-0.3)]<-1
mydata$votech2[mydata$votech/100>=(-0.3) & mydata$votech/100<(-0.2)]<-2
mydata$votech2[mydata$votech/100>=(-0.2) & mydata$votech/100<(-.1)]<-3
mydata$votech2[mydata$votech/100>=(-0.1) & mydata$votech/100<(-.05)]<-4
mydata$votech2[mydata$votech/100>=(-.05) & mydata$votech/100<0  ]<-5
mydata$votech2[is.na(mydata$votech/100)==TRUE]<-6
mydata$votech2[mydata$votech/100>=0 & mydata$votech/100<0.05]<-7
mydata$votech2[mydata$votech/100>=0.05 & mydata$votech/100<0.1]<-8
mydata$votech2[mydata$votech/100>=.1 & mydata$votech/100<.2]<-9
mydata$votech2[mydata$votech/100>=.2 & mydata$votech/100<.3]<-10
mydata$votech2[mydata$votech/100>=.3 ]<-11

mydata$votech2<-factor(mydata$votech2, ordered=TRUE)



p1 <- ggplot() +
  geom_polygon(data=mydata, aes(x=long, y=lat, group=group, fill=votech2), size=.001,colour="gray") +
  geom_path(data=states, aes(x=long, y=lat, group=group), color="black")+
  theme(axis.text=element_blank(),
              axis.ticks=element_blank(),
              axis.title.x=element_blank(),
              axis.title.y=element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.border = element_blank(),
              panel.background = element_blank())+
   					 theme(legend.position=c(.5,-.1),legend.text = element_text(size=7),legend.title = element_text(size=10),legend.key.size = unit(.3, "cm")) +
  scale_fill_brewer(name="Vote Change:", palette="RdBu",direction=-1,breaks=c(1,2,3,4,5,6,7,8,9,10,11),labels=c(">30% Dem",">20% Dem",">10% Dem",">5% Dem","<5% Dem","NA/Unopposed","<1% Rep","<5% Rep","<10% Rep","<20% Rep","30% Rep")) +guides(fill=guide_legend(ncol=3, title.position="top"))

p1+coord_map("albers",lat0=39, lat1=45)



ggsave("map_votech1.pdf", width=6.5, units="in")
```


Centering all variables:
```{r}

center<-function(myvar){
  (myvar-mean(myvar, na.rm=TRUE))
}

tr_cr<-trumpilot %>% 
  select(sales_USD12,prod_bush12, sales_USD17,prod_bush17,GDP2015,perc_HS_GED, perc_BA, totpop1317,urb2010,perclatino1317,percblack1317,perc_otherrace,perc_unemploy1317,perc_foreign1317,trumpmarg, Rep16perc,corn_bush,corn_usd,cotton_bales,cotton_usd)  %>% 
  mutate(ln_GDP2015=log(GDP2015/1000), ln_pop=log(totpop1317/1000), perc_latino1317=perclatino1317*100,perc_HS_GED=perc_HS_GED*100,perc_BA=perc_BA*100,perc_black1317=percblack1317*100,ln_perc_black1317=log(percblack1317*100+1),ln_perc_latino1317=log(perclatino1317*100+1),perc_unemploy1317=perc_unemploy1317*100,perc_foreign1317=perc_foreign1317*100,ln_perc_foreign1317=log(perc_foreign1317*100+1),sales_USD12=sales_USD12/100000, prod_bush12=prod_bush12/1000000,prod_bush17=prod_bush17/1000000, ln_sales_USD12=log(sales_USD12/1000+1), ln_bush_pop=log(prod_bush12/(1000000*totpop1317)+1), ln_gdpc=log(GDP2015/totpop1317), ln_prod_bush12=log(prod_bush12+1),ln_prod_bush17=log(prod_bush17+1),sales_USD17=sales_USD17/100000,ln_sales_USD17=log(sales_USD17/1000+1), Rep16perc=Rep16perc*100,trumpmarg=trumpmarg*100, corn_usd=as.numeric(as.character(corn_usd))) %>% mutate(ln_perc_other=log(100*perc_otherrace+1),perc_other=100*perc_otherrace,corn_bush=corn_bush/1000000,cotton_bales=cotton_bales/10000,ln_corn_bush=log(corn_bush+1), ln_corn_usd=log(corn_usd/1000+1),ln_cotton_bales=log(cotton_bales+1),ln_cotton_usd=log(cotton_usd/1000+1)) 

tr_cr<-do.call(cbind,lapply(tr_cr,center))
newnames<-paste(names(as.data.frame(tr_cr)),"_cr",sep="")
colnames(tr_cr)<-newnames


trumpilot1<-cbind.data.frame(subset(trumpilot, select=c("FIPS","statefips","state","county","Dist","Splitflag","votech","soycounty_bin")),tr_cr)
trumpilot1$ln_gdpc_sq_cr=trumpilot1$ln_gdpc_cr^2


```

Centering variables in non-split county subset:
```{r}

center<-function(myvar){
  (myvar-mean(myvar, na.rm=TRUE))
}

tr_cr<-trumpilot %>% filter(Splitflag==0) %>%
  select(sales_USD12,prod_bush12,sales_USD17,prod_bush17, GDP2015,perc_HS_GED, perc_BA, totpop1317,urb2010,perclatino1317,percblack1317,perc_unemploy1317,perc_foreign1317,trumpmarg, PVI2015,perc_otherrace, Rep16perc,corn_bush,corn_usd,cotton_bales,cotton_usd)  %>% 
  mutate(ln_GDP2015=log(GDP2015/1000),ln_pop=log(totpop1317/1000), perc_latino1317=perclatino1317*100,perc_HS_GED=perc_HS_GED*100,perc_BA=perc_BA*100,perc_black1317=percblack1317*100,ln_perc_black1317=log(percblack1317*100+1),ln_perc_latino1317=log(perclatino1317*100+1),perc_unemploy1317=perc_unemploy1317*100,perc_foreign1317=perc_foreign1317*100,ln_perc_foreign1317=log(perc_foreign1317*100+1),sales_USD12=sales_USD12/100000, prod_bush12=prod_bush12/1000000, ln_sales_USD12=log(sales_USD12/1000+1),prod_bush17=prod_bush17/1000000, ln_bush_pop=log(prod_bush12/(1000000*totpop1317)+1), ln_prod_bush17=log(prod_bush17+1),sales_USD17=sales_USD17/100000,ln_sales_USD17=log(sales_USD17/1000+1), ln_gdpc=log(GDP2015/totpop1317), ln_prod_bush12=log(prod_bush12+1),Rep16perc=Rep16perc*100,trumpmarg=trumpmarg*100,corn_usd=as.numeric(as.character(corn_usd))) %>% mutate(ln_perc_other=log(100*perc_otherrace+1),perc_other=100*perc_otherrace, corn_bush=corn_bush/1000000,cotton_bales=cotton_bales/10000,ln_corn_bush=log(corn_bush+1), ln_corn_usd=log(corn_usd/1000+1),ln_cotton_bales=log(cotton_bales+1),ln_cotton_usd=log(cotton_usd/1000+1)) 

tr_cr<-do.call(cbind,lapply(tr_cr,center))
newnames<-paste(names(as.data.frame(tr_cr)),"_cr",sep="")
colnames(tr_cr)<-newnames


trumpilot2<-cbind.data.frame(subset(trumpilot[trumpilot$Splitflag==0,], select=c("FIPS","statefips","state","county","Dist","votech","incshiftr","distname", "soycounty_bin")),tr_cr)
trumpilot2$ln_gdpc_sq_cr=trumpilot2$ln_gdpc_cr^2
```



Table 1:
```{r}
library(lme4)
#Model 1:
summary(m1<-lmer( votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr +  Rep16perc_cr + (1  | statefips),data = trumpilot1))

#Model 2:
summary(m2<-lmer(data = trumpilot1, votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr+ (1  | statefips)))

#Model 3
summary(m3<-lmer(data = trumpilot2, votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+  Rep16perc_cr + incshiftr + PVI2015_cr+ (1 | statefips)))

#Model 4
summary(m4<-lmer(data = trumpilot2, votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+ Rep16perc_cr + incshiftr + PVI2015_cr+(1 | statefips)))


```


Online Appendix:
```{r}
#Figure 2: 2012 Soy Production (bushels)
mydata$soy_scale1<-0
mydata$soy_scale1[mydata$prod_bush12>=62 & mydata$prod_bush12<105242]<-1
mydata$soy_scale1[mydata$prod_bush12>=105242 & mydata$prod_bush12<703941]<-2
mydata$soy_scale1[mydata$prod_bush12>=703941 & mydata$prod_bush12<2339635]<-3
mydata$soy_scale1[mydata$prod_bush12>=2339635 ]<-4
mydata$soy_scale1<-factor(mydata$soy_scale1, ordered=TRUE)

p <- ggplot() +
  geom_polygon(data=mydata, aes(x=long, y=lat, group=group, fill=soy_scale1),color="gray",  size=.0001) +
  geom_path(data=states, aes(x=long, y=lat, group=group), color="black")+
  theme(axis.text=element_blank(),
              axis.ticks=element_blank(),
              axis.title.x=element_blank(),
              axis.title.y=element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.border = element_blank(),
              panel.background = element_blank())+
   					 theme(legend.position=c(.5,-0.05),legend.text = element_text(size=7),legend.title = element_text(size=10),legend.key.size = unit(.3, "cm")) +
  scale_fill_brewer(name=" ", palette="Reds", labels=c("none","<25 pctl","< median","< 75 pctl","> 75 pctl")) +guides(fill=guide_legend(nrow=2, title.position="top"))

p +coord_map("albers",lat0=39, lat1=45)

ggsave("soybean_map2.pdf", width=6.5, units="in")
```

Figure 3: Expected Change in Vote Share of the Two Major Parties as a Function of Soy Production
```{r}
library(MASS)
theme_set(theme_grey() + theme(panel.background = element_rect(fill = NA, color = 'black'))+ theme(axis.text=element_text(size=10),
					axis.title=element_text(size=12,face="bold")))

ln_prod_bush12=seq(min(log(trumpilot$prod_bush12+1), na.rm=TRUE),max(log(trumpilot$prod_bush12+1), na.rm=TRUE), by=.1)
mycoef<-as.matrix(fixef(m1)[1:2])
mycov<-vcov(m1)



#A rural non-diverse county:
mycoef<-as.matrix(fixef(m1))
urb2010_cr<-min(trumpilot1$urb2010_cr, na.rm=TRUE)
ln_gdpc_cr<-mean(trumpilot1$ln_gdpc_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
ln_gdpc_sq_cr<-ln_gdpc_cr^2
perc_unemploy1317_cr<-mean(trumpilot1$perc_unemploy1317_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
perc_HS_GED_cr<-mean(trumpilot1$perc_HS_GED_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
perc_BA_cr<-mean(trumpilot1$perc_BA_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
percblack1317_cr<-mean(trumpilot1$ln_perc_black1317_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
perclatino1317_cr<-mean(trumpilot1$ln_perc_latino1317_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
perc_other_cr<-mean(trumpilot1$ln_perc_other_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
perc_foreign1317_cr<-mean(trumpilot1$ln_perc_foreign1317_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)
Rep16perc_cr<-mean(trumpilot1$Rep16perc_cr[trumpilot1$urb2010_cr==min(trumpilot1$urb2010_cr, na.rm=TRUE)], na.rm=TRUE)

yhat<-cbind(1,ln_prod_bush12,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,percblack1317_cr,perc_other_cr,perclatino1317_cr,perc_foreign1317_cr,Rep16perc_cr)%*%mycoef

mycoef1000<-mvrnorm(1000,mycoef,mycov)

yhat1000<-function(mycoef1000){
yhat1000<-cbind(1,ln_prod_bush12,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,percblack1317_cr,perc_other_cr,perclatino1317_cr,perc_foreign1317_cr,Rep16perc_cr)%*%mycoef1000
}

sims<-yhat1000(mycoef1000[1,])
sims1000<-apply(mycoef1000,1, function(x) yhat1000(x))
yhat_sim<-rowMeans(sims1000)
se_sim<-apply(sims1000, 1, sd)

p<-ggplot()+geom_line(aes(x=(exp(ln_prod_bush12)/1000),y=yhat_sim))+geom_ribbon(aes(x=(exp(ln_prod_bush12)/1000), ymin=(yhat_sim-2*se_sim),ymax=(yhat_sim+2*se_sim)), alpha=.1) 
p+ geom_rug(aes(x=trumpilot$prod_bush12/1000),alpha = .5)


p+ geom_jitter(data=subset(trumpilot,subset=!is.na(votech)),aes(x=(prod_bush12/1000+0.001),y=votech),alpha=0)+scale_y_continuous("Expected Change in Republican Vote Share, %",limits=c(-50,0),breaks=seq(from=-50,to=0,by=10), labels=seq(from=-50,to=0,10))+scale_x_log10("Soy Production, Thousands of Bushels", breaks=c(0.001,1,10,100,1000,10000), labels=c("none","1","10","100","1000","10000"))+geom_rug(data=subset(trumpilot,subset=!is.na(votech)),aes(x=(prod_bush12/1000+0.001),y=votech),alpha = .2, position  = position_jitter(width = 0.1, height = 0.1)) 
p 
ggsave("fitted2.pdf")


```

Table 2: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function
of Soy Production, with State-Level Fixed Effects
```{r}
detach("package:lme4", unload=TRUE) #MASS in lme4 clashes with dlyr's select function
detach("package:MASS", unload=TRUE)

trumpilot1<-left_join(trumpilot1,fixeffs, by="FIPS")
trumpilot2<-left_join(trumpilot2,fixeffs, by="FIPS")

#Model 1:
myvars_stateeffs<- trumpilot1 %>% select(votech,ln_prod_bush12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,statefips)
summary(m1_sfe<-lm(votech ~ .,data = myvars_stateeffs))

#Model 2:
myvars_stateeffs2<- trumpilot1 %>% select(votech,ln_sales_USD12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,statefips)
summary(m2_sfe<-lm(votech ~ .,data = myvars_stateeffs2))

#Model 3:
myvars_stateeffs3<- trumpilot2 %>% select(votech,ln_prod_bush12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,incshiftr,PVI2015_cr,statefips)
summary(m3_sfe<-lm(votech ~ .,data = myvars_stateeffs3))

#Model 4:
myvars_stateeffs4<- trumpilot2 %>% select(votech,ln_sales_USD12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr, incshiftr,PVI2015_cr,statefips)
summary(m4_sfe<-lm(votech ~ .,data = myvars_stateeffs4))

```

Table 3: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, with District-Level Fixed Effects
```{r}

#M1
myvars_disteffs<- trumpilot1 %>% select(votech,ln_prod_bush12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,fixeff_AL2:fixeff_WY1)
summary(m1_dfe<-lm(votech ~ .,data = myvars_disteffs))

#M2
myvars_disteffs2<- trumpilot1 %>% select(votech,ln_sales_USD12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,fixeff_AL2:fixeff_WY1)
summary(m2_dfe<-lm(votech ~ .,data = myvars_disteffs2))

#M3
myvars_disteffs3<- trumpilot2 %>% select(votech,ln_prod_bush12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,incshiftr,PVI2015_cr,fixeff_AL2:fixeff_WY1)
summary(m3_dfe<-lm(votech ~ .,data = myvars_disteffs3))

#M4
myvars_disteffs4<- trumpilot2 %>% select(votech,ln_sales_USD12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,incshiftr,PVI2015_cr,fixeff_AL2:fixeff_WY1)
summary(m4_dfe<-lm(votech ~ .,data = myvars_disteffs4))



```

Table 4: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, with State- and District-Level Fixed Effects
```{r}
#M1
myvars_both<- trumpilot1 %>% select(votech,ln_prod_bush12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,statefips,fixeff_AL2:fixeff_WY1)
summary(m1_2fe<-lm(votech ~ .,data = myvars_both))

#M2
myvars_both2<- trumpilot1 %>% select(votech,ln_sales_USD12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,statefips,fixeff_AL2:fixeff_WY1)
summary(m2_2fe<-lm(votech ~ .,data = myvars_both2))


#M3
myvars_both3<- trumpilot2 %>% select(votech,ln_prod_bush12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,incshiftr,PVI2015_cr,statefips,fixeff_AL2:fixeff_WY1)
summary(m3_2fe<-lm(votech ~ .,data = myvars_both3))

#M4
myvars_both4<- trumpilot2 %>% select(votech,ln_sales_USD12_cr,ln_gdpc_cr,ln_gdpc_sq_cr,perc_unemploy1317_cr,perc_HS_GED_cr,perc_BA_cr,urb2010_cr,ln_perc_black1317_cr,ln_perc_other_cr,ln_perc_latino1317_cr,ln_perc_foreign1317_cr,Rep16perc_cr,statefips,incshiftr,PVI2015_cr,fixeff_AL2:fixeff_WY1)
summary(m4_2fe<-lm(votech ~ .,data = myvars_both4))
```

Table 5: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Accounting for Corn and Cotton Production
```{r}
library(lme4)
#M1
summary(m1_othercrops<-lmer( votech ~ ln_prod_bush12_cr + ln_corn_bush_cr+ln_cotton_bales_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + (1 | statefips),data = trumpilot1))

#M2
summary(m2_othercrops<-lmer(data = trumpilot1, votech ~ ln_sales_USD12_cr +ln_corn_usd_cr+ln_cotton_usd_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr+ (1  | statefips)))

#M3
summary(m3_othercrops<-lmer(data = trumpilot2, votech ~ ln_prod_bush12_cr + ln_corn_bush_cr+ln_cotton_bales_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + incshiftr + PVI2015_cr+ (1 | statefips)))

#M4
summary(m4_othercrops<-lmer(data = trumpilot2, votech ~ ln_sales_USD12_cr +ln_corn_usd_cr+ln_cotton_usd_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + incshiftr + PVI2015_cr  + (1 | statefips)))

```

Table 6: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Controlling for Trump’s 2016 Margin

```{r}
#M1
summary(m1_trump<-lmer( votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+trumpmarg_cr + (1  | statefips),data = trumpilot1))

#M2
summary(m2_trump<-lmer(data = trumpilot1, votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+trumpmarg_cr+ (1  | statefips)))

#M3
summary(m3_trump<-lmer(data = trumpilot2, votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+trumpmarg_cr + incshiftr + PVI2015_cr+ (1 | statefips)))

#M4
summary(m4_trump<-lmer(data = trumpilot2, votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+trumpmarg_cr + incshiftr + PVI2015_cr+(1 | statefips)))
```

Table 7: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Accounting for Spatial Dependence in Soy Production

```{r}
library(spdep)
library(spatialreg)
W<-read.table("contigmat.txt", row.names=1, header=TRUE)
#Row-standardize:
W_st<-W
W_st[rowSums(W)!=0]<-W[rowSums(W)!=0]/rowSums(W)[rowSums(W)!=0]
trumpilot1$spatX<-t(trumpilot1$ln_prod_bush12_cr%*%as.matrix(W_st))
trumpilot1$spatX1<-t(trumpilot1$ln_sales_USD12_cr%*%as.matrix(W_st))
trumpilot2<-left_join(trumpilot2,subset(trumpilot1, select=c("FIPS","spatX","spatX1")),by=c("FIPS"))

#Spatial X model:
#M1
summary(m1_spX<-lmer( votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr  +spatX+ (1  | statefips),data = trumpilot1))


#M2
summary(m2_spX<-lmer(data = trumpilot1, votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+  Rep16perc_cr+ spatX1+ (1  | statefips)))

#M3
summary(m3_spX<-lmer(data = trumpilot2, votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+  Rep16perc_cr+ spatX + incshiftr + PVI2015_cr+ (1 | statefips)))

#M4
summary(m4_spX<-lmer(data = trumpilot2, votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+  Rep16perc_cr + spatX1+ incshiftr + PVI2015_cr+(1 | statefips)))


```

Table 8: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Accounting for Spatial Dependence in Vote
```{r}
W1<-mat2listw(as.matrix(W), row.names=as.integer(as.character(trumpilot1$FIPS)))
W2<-nb2listw(W1$neighbours, glist=NULL, style="W", zero.policy=TRUE)

#M1:
spY1<-lagsarlm(votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr +statefips , data=trumpilot1, W2, zero.policy=TRUE)
summary(spY1)


#M2
spY2<-lagsarlm(votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr +statefips , data=trumpilot1, W2, zero.policy=TRUE)
summary(spY2)

```

Table 9: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Accounting for Spatial Dependence in the Error Term
```{r}
#M1
spE1<-errorsarlm(votech ~ ln_prod_bush12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr +statefips, data=trumpilot1, W2, zero.policy=TRUE)
summary(spE1)


#M2
spE2<-errorsarlm(votech ~ ln_sales_USD12_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr +statefips , data=trumpilot1, W2, zero.policy=TRUE)
summary(spE2)
```

Matching:
```{r}
library(cem)
library(optmatch)
library(RItools)

trumpilot$soycounty_bin<-as.numeric(trumpilot$prod_bush12>10000)
trumpilot4<-trumpilot %>% select(FIPS,votech,soycounty_bin,GDP2015,totpop1317,perc_unemploy1317,perc_HS_GED,perc_BA,urb2010,perclatino1317,percblack1317,perc_otherrace,perc_foreign1317,Rep16perc,votech,prod_bush12,sales_USD12,statefips) %>% 
mutate(ln_gdpc=log(GDP2015/totpop1317),                                                         perc_latino1317=perclatino1317*100,ln_perc_latino1317=log(perclatino1317*100+1),                             perc_HS_GED=perc_HS_GED*100,                                                                                 perc_BA=perc_BA*100,                                               ln_perc_black1317=log(percblack1317*100+1),
perc_unemploy1317=perc_unemploy1317*100,                                                                     perc_foreign1317=perc_foreign1317*100,ln_perc_foreign1317=log(perc_foreign1317*100+1),        sales_USD12=sales_USD12/100000, ln_sales_USD12=log(sales_USD12/1000+1),                                      prod_bush12=prod_bush12/1000000,ln_prod_bush12=log(prod_bush12+1),                                           Rep16perc=Rep16perc*100,
ln_perc_other=log(100*perc_otherrace+1)) %>%
select(statefips,soycounty_bin,votech,ln_gdpc,ln_perc_latino1317,perc_HS_GED,perc_BA,urb2010,ln_perc_black1317,perc_unemploy1317,ln_perc_foreign1317,ln_sales_USD12,ln_prod_bush12,Rep16perc,ln_perc_other,FIPS) 
row.names(trumpilot4)<-trumpilot4$FIPS
trumpilot4<-na.omit(trumpilot4)
trumpilot4<-trumpilot4 %>% select(statefips,soycounty_bin,votech,ln_prod_bush12,ln_gdpc,ln_perc_latino1317,perc_HS_GED,perc_BA,urb2010,ln_perc_black1317,perc_unemploy1317,ln_perc_foreign1317,Rep16perc,ln_perc_other) 

#Balance in the original data:
xBalance(soycounty_bin ~  ln_gdpc+perc_unemploy1317+perc_HS_GED+perc_BA+urb2010+ln_perc_black1317+ln_perc_other+ln_perc_latino1317+ln_perc_foreign1317+Rep16perc, 
  data = trumpilot4, 
  report = c("chisquare.test","std.diffs"))


unemploycuts<-seq(from=24,to=87,by=7)
repcuts<-c(25,50,75)
gdpcuts<-quantile(trumpilot4$ln_gdpc, probs=c(.1,.2,.3,.4,.5,.6,.7,.8,.9))
HScuts<-c(10,30,40,50)
BAcuts<-c(10,15,30)
foreigncuts<-c(1000,5000)
trumpilot4$urbcut=as.factor(as.numeric(trumpilot4$urb2010>0))
mat1 <- cem(treatment = "soycounty_bin",  drop=c("votech","ln_prod_bush12","statefips","urb2010","ln_gdpc_sq"), cutpoints=list(ln_gdpc=gdpcuts,perc_unemploy1317=unemploycuts,perc_BA=BAcuts,ln_perc_black1317=3,
ln_perc_foreign1317=foreigncuts,ln_perc_other=log(11),Rep16perc=repcuts),data = trumpilot4, keep.all = TRUE)
mat1
mat1$breaks

#Table 10: Balance in the Post-Matched Sample, Measured as the Difference in Means Between the Treated and Control Groups
xBalance(soycounty_bin ~  ln_gdpc+perc_unemploy1317+perc_HS_GED+perc_BA+urb2010+ln_perc_black1317+ln_perc_other+ln_perc_latino1317+ln_perc_foreign1317+Rep16perc, 
  data = trumpilot4[mat1$matched==TRUE,], 
  report = c("chisquare.test","adj.mean.diffs","std.diffs"))


```


Figure 4: Balance in the Post-Matched Sample, Measured as the Overlap in the Empirical Density Distributions of the Treatment and Control Groups
```{r}
#Plot distributions:
data_long <- trumpilot4[mat1$matched==TRUE,] %>%
gather(covar, value, c("ln_gdpc","perc_unemploy1317","perc_HS_GED","perc_BA","urb2010","ln_perc_black1317","ln_perc_other","ln_perc_latino1317","ln_perc_foreign1317","Rep16perc"))
labels <- c(ln_gdpc="GDP/cap",perc_unemploy1317="Unemployed",perc_HS_GED="HS/GED",perc_BA="BA",urb2010="Urb",ln_perc_black1317="Black",ln_perc_other="Other Races",ln_perc_latino1317="Latino",ln_perc_foreign1317="Foreign",Rep16perc="Rep 16")

p<-ggplot(data=data_long)+geom_density(aes(x=value,linetype=as.factor(soycounty_bin)))+
facet_wrap(~covar, nrow=2,scales = "free",labeller=labeller(covar = labels))+
  theme(legend.position="none",strip.text = element_text(size=8),
          strip.background = element_rect( fill="white"))+ scale_y_continuous(breaks=NULL) + ylab("")+ scale_x_continuous(breaks=NULL) + xlab("")+scale_linetype_manual(values=c("solid","dashed"))
p

ggsave("matchplot.jpg")


```


Centering  variables in the matched sample:
```{r}

center<-function(myvar){
  (myvar-mean(myvar, na.rm=TRUE))
}

tr_cr<-trumpilot4[mat1$matched==TRUE,] %>% 
  select(ln_prod_bush12,ln_gdpc,ln_perc_latino1317,perc_HS_GED,perc_BA,urb2010,ln_perc_black1317,perc_unemploy1317,ln_perc_foreign1317,Rep16perc,ln_perc_other)  

tr_cr<-do.call(cbind,lapply(tr_cr,center))
newnames<-paste(names(as.data.frame(tr_cr)),"_cr",sep="")
colnames(tr_cr)<-newnames


trumpilot5<-cbind.data.frame(subset(trumpilot4[mat1$matched==TRUE,], select=c("statefips","votech","soycounty_bin")),tr_cr)
trumpilot5$ln_gdpc_sq_cr=trumpilot5$ln_gdpc_cr^2
trumpilot5$ln_gdpc_cu_cr=trumpilot5$ln_gdpc_cr^3

```

Table 11: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Estimated on a Matched Sample
```{r}
summary(mmatch1<-lmer(votech ~ ln_prod_bush12_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + (1 |statefips), 
  data = trumpilot5))
```

Table 12: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Estimated Using the Difference in Soy Production
```{r}
trumpilot1$ln_soy_bush_diff=log(trumpilot1$prod_bush17_cr-trumpilot1$prod_bush12_cr+1.5)
trumpilot2$ln_soy_bush_diff=log(trumpilot2$prod_bush17_cr-trumpilot2$prod_bush12_cr+1.67)



#M1
summary(m1_dd<-lmer(votech ~ ln_soy_bush_diff+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + (1 | statefips),data = trumpilot1[trumpilot1$soycounty_bin==1,]))


#M2
summary(m3_dd<-lmer(data = trumpilot2[trumpilot2$soycounty_bin==1,], votech ~ ln_soy_bush_diff + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr  + incshiftr + PVI2015_cr + (1 | statefips)))


```

Table 13: Change in Republican Vote Share Between 2016 and 2018 Elections as a Function of Soy Production, Using 2017 Data on Soy Production

```{r}
#M1
summary(m1_17<-lmer(votech ~ ln_prod_bush17_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + (1 | statefips),data = trumpilot1))


#M2
summary(m2_17<-lmer(votech ~ ln_sales_USD17_cr+ ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + (1 | statefips),data = trumpilot1))


#M3
summary(m3_17<-lmer(data = trumpilot2, votech ~ ln_prod_bush17_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr  + incshiftr + PVI2015_cr + (1 | statefips)))

#M4
summary(m4_17<-lmer(data = trumpilot2, votech ~ ln_sales_USD17_cr + ln_gdpc_cr+ln_gdpc_sq_cr+perc_unemploy1317_cr+perc_HS_GED_cr+perc_BA_cr+urb2010_cr+ln_perc_black1317_cr+ln_perc_other_cr+ln_perc_latino1317_cr+ln_perc_foreign1317_cr+Rep16perc_cr + incshiftr + PVI2015_cr + (1 | statefips)))
```

